Introduction à la modélisation statistique bayésienne

Thierry Phénix & Ladislas Nalborczyk
21 mai 2019

Planning

Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton

Plan

  • Rappels et objectifs
  • Introduction
  • Evaluation de la distribution posterieure par échantillonnage
  • Vérification du résultat

On a vu quoi jusqu'ici?

  • The first idea is that Bayesian inference is reallocation of credibility across possibilities.
  • The second foundational idea is that the possibilities, over which we allocate credibility, are parameter values in meaningful mathematical models.
    Kruschke (2015)


3 configurations possibles:

  • prior et likelihood sont conjugées l'une de l'autre -> expression analytique de la distribution postérieure
  • méthode grille -> Calcul de la valeur exacte de la postérieure en un nombre fini de points
  • échantillonner un grand nombre de points dans l'espace des paramètres

Les objectifs du cours


\( \longrightarrow~ \) Présenter le principe de base de l'échantillonnage : les MCMC

\( \longrightarrow~ \) Présenter les méthodes classiques (point de vue historique)

\( \longrightarrow~ \) Montrer les forces mais aussi les faiblesses de ces méthodes

\( \longrightarrow~ \) Donner des outils de contrôle sur ces méthodes

\( \longrightarrow~ \) Appliquer ces méthodes à un cas simple

Planning

Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton

Plan

  • Rappels et objectifs
  • Introduction
  • Evaluation de la distribution posterieure par échantillonnage
  • Vérification du résultat

MCMC, qu'y-a-t-il dans ce nom?

  • Markov chain Monte Carlo
    \( \longrightarrow~ \) Échantillonnage aléatoire
    \( \longrightarrow~ \) Le résultat est un ensemble de valeurs du paramètres

  • Markov chain Monte Carlo
    \( \longrightarrow~ \) Les valeurs sont générées sous forme de séquence (liaison de dépendance)
    \( \longrightarrow~ \) Indice temporel pour identifier la place dans la chaine
    \( \longrightarrow~ \) Le résultat est de la forme \( ~:~\theta^1, \theta^2,\theta^3, \cdots, \theta^t \)

  • Markov chain Monte Carlo
    \( \longrightarrow~ \) La valeur de paramètre générée ne dépend que de la valeur du paramètre précédent \( \longrightarrow~P(\theta^{t+1}~|~\theta^{t},\theta^{t-1}, \cdots,\theta^{1})~=~P(\theta^{t+1}~|~\theta^{t}) \)

Méthodes Monte Carlo

Le terme méthode de Monte-Carlo, désigne une famille d'algorithmiques visant à calculer une valeur numérique approchée en utilisant des procédés aléatoires, c'est-à-dire des techniques probabilistes.
Wikipedia

Cette méthode a été inventé en 1947 par Nicholas Metropolis, et publié pour la première fois en 1949 dans un article coécrit avec Stanislaw Ulam.

Jacques Bernoulli
Nicolas Metropolis (1915-1999)

Méthodes Monte Carlo

Les méthodes Monte Carlo designe une famille d'algorithmes qui ont pour but de calculer des valeurs numériques approchées à partir de procédés aléatoires


\( \longrightarrow~ \) Utiliser une méthode Monte Carlo pour calculer une valeur approchée de \( \pi \)
La probabilité que le point \( M \) appartienne au disque est \( \pi⁄4 \).

Jacques Bernoulli

https://en.wikipedia.org/wiki/Monte_Carlo_method

Méthodes Monte Carlo

Les méthodes Monte Carlo designe une famille d'algorithmes qui ont pour but de calculer des valeurs numériques approchées à partir de procédés aléatoires


\( \longrightarrow~ \) Utiliser une méthode Monte Carlo pour chercher un maximum (Optimisation)
Par exemple déterminer le maximum d'une fonction.

Algorithme Simulated annealing

https://en.wikipedia.org/wiki/Simulated_annealing

Méthodes Monte Carlo

Les méthodes Monte Carlo designe une famille d'algorithmes qui ont pour but de calculer des valeurs numériques approchées à partir de procédés aléatoires


\( \longrightarrow~ \) Utiliser une méthode Monte Carlo pour approximer une distribution de probabilité

On connait les priors \( ~P(\theta_1) \) et \( ~P(\theta_2) \)
On connait la likelihood \( ~P(D~|~\theta_1, \theta_2) \)
Mais on ne sait pas calculer la distribution postéreure \( ~P(\theta_1, \theta_2~|~D)=\frac{P(D~|~\theta_1, \theta_2)~P(\theta_1)P(\theta_1)~}{P(D)} \)

DONC on sait calculer la distribution postérieure à une constante près car \( P(D) \) est une constante.
On va évaluer la distribution postérieure à partir de notre connaissance sur les priors et la likelihood avec une méthode d'échantillonnage

Méthodes Monte Carlo : Exemple

Considérons un exemple simple

Soit un paramètre \( \theta \) qui ne prend que 7 valeurs, et une répartition de ce paramètre :

Question : Quelle est la distribution de \( \theta~ \)?


Remarque : mais c'est évident… ICI

Estimation de cette distribution par tirage aléatoire

On tire aléatoirement un grand nombre de point au hasard parmi ces 28 cases!
(Comme pour le calcul de \( \pi \))

Le résultat du tirage \( (4,6,2,3,7,7,5,3,6,5,1, \dots) \)

Méthodes Monte Carlo : Exemple


plot of chunk unnamed-chunk-2


Que peut-on dire de cette méthode dans le cadre d'une distribution :
\( ~\bullet~ \) Elle converge vers la distribution
\( ~\bullet~ \) Aucun contrôle sur la vitesse de convergence
\( ~\bullet~ \) Nécessite beaucoup de points (7 valeurs de paramètres et 100 tirages)

\( ~\rightarrow~ \) En l'état, résultat pas concluant !

Planning


Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton


Plan

  • Rappels et objectifs
  • Introduction
  • Evaluation de la distribution posterieure par échantillonnage
    • L'algorithme Metropolis
    • L'algorithme Gibbs
    • L'algorithme Hamilton
  • Vérification du résultat

MCMC : algorithme Metropolis

Cet algorithme a été présenté pour la première fois en 1953
Les auteurs :
Nicholas Metropolis, Arianna W. Rosenbluth, Marshall Rosenbluth, Augusta H. Teller, et Edward Teller


Marshall Nicholas Rosenbluth
(5 February 1927 – 28 September 2003)
Physicien Americain (spécialiste du plasma)

MCMC : algorithme Metropolis

Première méthode
Le problème n'est pas la convergence, mais la vitesse à laquelle la méthode converge
Pour augmenter la vitesse de convergence, il faudrait faciliter l'accès aux valeurs de paramètres les plus représentées.

Principe
\( \bullet~ \) On fait une proposition de déplacement sur la base de la valeur courante du paramètre
\( \bullet~ \) on réalise un tirage aléatoire pour accepter ou rejeter la nouvelle position

2 Idées
\( \bullet~ \) la proposition doit favoriser les valeurs de paramètre les plus probables
\( ~~\longrightarrow~ \) On parcourt plus souvent ces valeurs de paramètres

\( \bullet~ \) la proposition doit se limiter aux valeurs adjacentes du paramètre courant
\( ~~\longrightarrow~ \) On augmente la vitesse de convergence en restant là où se trouve l'information

MCMC : algorithme Metropolis

Algorithme de Metropolis

\( \longrightarrow~ \) Sélectionner un point de départ
On peut sélectionner n'importe quelle valeur

MCMC : algorithme Metropolis

Algorithme de Metropolis

\( \longrightarrow~ \) Sélectionner un point de départ
On peut sélectionner n'importe quelle valeur

\( \longrightarrow~ \) Faire une proposition de déplacement
Faire un tirage et proposer un déplacement

MCMC : algorithme Metropolis

Algorithme de Metropolis

\( \longrightarrow~ \) Sélectionner un point de départ
On peut sélectionner n'importe quelle valeur

\( \longrightarrow~ \) Faire une proposition de déplacement
Faire un tirage et proposer un déplacement

\( \longrightarrow~ \) Accepter ou rejeter la proposition de déplacement
Faire un tirage suivant la probabilité :

\[ P_{déplacement}=min \left(\frac{P(\theta_{proposé})}{P(\theta_{courant})},1 \right) \]

Remarque : Plus la valeur de \( \theta \) est grande et plus petite est la probabilité de rester sur place

MCMC : algorithme Metropolis

Algorithme de Metropolis

\( \longrightarrow~ \) Sélectionner un point de départ
On peut sélectionner n'importe quelle valeur

\( \longrightarrow~ \) Faire une proposition de déplacement
Faire un tirage et proposer un déplacement

\( \longrightarrow~ \) Accepter ou rejeter la proposition de déplacement
Faire un tirage suivant la quasi probabilité :

\[ P_{déplacement}=min \left(\frac{P(\theta_{proposé})}{P(\theta_{courant})},1 \right) \]

\( \longrightarrow~ \) La position calculée devient la nouvelle position de départ

Méthodes Monte Carlo vs Algorithme Metropolis

plot of chunk unnamed-chunk-3

plot of chunk unnamed-chunk-4

MCMC : algorithme Metropolis

Application au lancer de pièce (cas continu)

\( ~\bullet~ \) La likelihood est donnée par : \( \color{orangered}{P(z~|~\theta,N)=\theta^z~(1-\theta)^{(N-z)}} \)
\( ~\bullet~ \) Le prior est donné par : \( \color{steelblue}{P(\theta~|~a,b) \propto \theta^{(a-1)}(1-\theta)^{(b-1)}} \)
\( ~\bullet~ \) Le paramètre que l'on cherche à estimer \( \theta \) prend ses valeurs dans l'intervalle \( \left[0,1\right] \)


Problème 1: Comment définir la proposition de déplacement?


Le déplacement est modélisé par une distribution normale :\( ~\Delta\theta~\sim ~\mathcal{N}0,\sigma) \)
\( \longrightarrow~ \) la moyenne \( \mu \) vaut \( 0 \) : le déplacement se fait autour de la valeur courante du paramètre
\( \longrightarrow~ \) La variance reste à déterminer, elle contrôle l'éloignement de la nouvelle valeur

MCMC : algorithme Metropolis

Problème 2: Quelle probabilité utiliser pour accepter ou refuser le déplacement?

Nous utilisons le produit de la likelihood et du prior : \( P(\theta)~\propto~ \color{orangered}{\theta^z~(1-\theta)^{(N-z)}}\color{steelblue}{\theta^{(a-1)}(1-\theta)^{(b-1)}} \)

La probabilité d'accepter le déplacement : \( ~P_{déplacement}=min \left(\frac{P(\theta_{cur}+\Delta\theta)}{P(\theta_{cur})},1 \right) \)

REMARQUE : le rapport \( \frac{P(\theta_{cur}+\Delta\theta)}{P(\theta_{cur})} \) est le même que l'on utilise la distribution postérieure ou le produit prior par likelihood !

MCMC : algorithme Metropolis

L'algorithme Metropolis

\( \longrightarrow~ \) Sélectionner un point de départ
\( ~\bullet~ \) Il faut choisir \( ~\theta \in \left[0,1\right] \)
\( ~\bullet~ \) Seule contrainte \( ~P(\theta_{init}) \ne 0 \)

\( \longrightarrow~ \) Faire une proposition de déplacement
\( ~\bullet~ \) Faire un tirage suivant \( ~\mathcal{N}0,\sigma) \)

\( \longrightarrow~ \) Accepter ou rejeter la proposition de déplacement
Faire un tirage suivant la probabilité :

\[ P_{déplacement}=min \left(\frac{P(\theta_{cur}+\Delta\theta)}{P(\theta_{cur})},1 \right) \]

\( \longrightarrow~ \) La position calculée devient la nouvelle position de départ

plot of chunk unnamed-chunk-5

MCMC : algorithme Metropolis

Le choix de sigma dans la proposition de déplacement

Deux indice permettre d'évaluer la qualité de l'échantillonnage :

\( \rightarrow~ \) Le rapport entre le nombre de déplacements proposés et le nombre de déplacements accéptés

\( \rightarrow~ \) L'effective sample size
C'est le nombre de déplacement qui ne sont pas corrélés avec les précédents

MCMC : algorithme Metropolis

Le choix de sigma dans la proposition de déplacement



\( \rightarrow~ \) Toutes les propositions de déplacement (ou presque) sont acceptées

\( \rightarrow~ \) Peu de valeurs effectives

Il va falloir beaucoup d'itérations pour avoir un résultat satisfaisant

MCMC : algorithme Metropolis

Le choix de sigma dans la proposition de déplacement



\( \rightarrow~ \) Les propositions de déplacement sont rarement acceptées

\( \rightarrow~ \) Peu de valeurs effectives

Il va falloir beaucoup d'itérations pour avoir un résultat satisfaisant

Planning


Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton


Plan

  • Rappels et objectifs
  • Introduction
  • Evaluation de la distribution posterieure par échantillonnage
    • L'algorithme Metropolis
    • L'algorithme Gibbs
    • L'algorithme Hamilton
  • Vérification du résultat

MCMC : algorithme Gibbs

Cet algorithme a été présenté pour la première fois en 1984
par les frêres Geman


https://en.wikipedia.org/wiki/Donald_Geman

https://en.wikipedia.org/wiki/Stuart_Geman

MCMC : algorithme Gibbs - modèle à deux variables

L'algorithme Gibbs pour évaluer une distribution posterieure ne fonctionne qu'avec des modèles composés d'au moins deux variables.


Définition du modèle bivarié

On considère une situation dans laquelle on dispose de deux pièces de monaie.
A chaque pièce on associe un paramètre ou biais \( ~\theta_1 \) et \( ~\theta_2 \) (la probabilité d'obtenir \( Face \)).
On cherche à estimer les biais de ces deux pièces : \( ~P(\theta_1, \theta_2~|~D ) \)
Nos variables \( \theta_1 \) et \( \theta_2 \) sont indépendantes et identiquement distribuées (iid)

MCMC : modèle à deux variables

Définition du piror

On considère que chacun des biais est décrit par une distribution beta :
\( \theta_1 \sim beta(\theta_1~|~a_1,b_1) \) et \( ~\theta_2 \sim beta(\theta_2~|~a_2,b_2) \)

On peut donc calculer la distribution a priori de \( ~P(\theta_1, \theta_2) \) :

\[ \begin{align} P(\theta_1, \theta_2) &= P(\theta_1)~P(\theta_2) \\ &= \theta_1^{(a_1-1)}~(1-\theta_1)^{(b_1-1)}\theta_2^{(a_2-1)}~(1-\theta_2)^{(b_2-1)} \end{align} \]


Représentation graphique de la distribution

MCMC : modèle à deux variables

Description des données

On dispose d'un ensemble de lancers pour les deux pièces
On appelle \( ~D_1 \) les données issues de la première pièce, composé de \( ~z_1 \) Face sur \( ~N_1 \) lancers.
On appelle \( ~D_2 \) les données issues de la seconde pièce, composé de \( ~z_2 \) Face sur \( ~N_2 \) lancers.

On note les données : \( D = \{ z_1, N_1, z_2, N_2 \} \)


Définition de la likelihood


\[ \begin{align} P(D~|~\theta_1,\theta_2) &= \prod_{y_{1i} \in D_1} P(y_{1i}~|~\color{orangered}{\theta_1})~ \prod_{y_{2i} \in D_2} P(y_{2i}~|~\color{orangered}{\theta_2}) &&\mbox{Règle d'indépendance}\\ &= \theta_1^{z_1}~(1-\theta_1)^{(N1-z_1)}\theta_2^{z_2}~(1-\theta_2)^{(N2-z_2)}~ \end{align} \]

MCMC : modèle à deux variables

Résolution analytique


\[ \begin{align} \color{purple}{P(\theta_1,\theta_2~|~D)} &= \color{orangered}{P(D~|~\theta_1,\theta_2)}~\color{steelBlue}{P(\theta_1,\theta_2)}/\color{green}{P(D)} \\ &= \color{orangered}{\theta_1^{z_1}~(1-\theta_1)^{(N1-z_1)}\theta_2^{z_2}~(1-\theta_2)^{(N2-z_2)}}~\color{steelBlue}{P(\theta_1,\theta_2)}/\color{green}{P(D)} \\ &= \frac{\theta_1^{(\color{orangered}{z_1}+\color{steelBlue}{a_1-1})}(1-\theta_1)^{(\color{orangered}{N1-z_1}+\color{steelBlue}{b_1-1})}\theta_2^{(\color{orangered}{z_2}+\color{steelBlue}{a_2-1})}(1-\theta_2)^{(\color{orangered}{N2-z_2}+\color{steelBlue}{b_2-1})}~}{\color{green}{P(D)}~\color{steelBlue}{B(a_1,b_1)~B(a_2,b_2)}} \end{align} \]

Donc finalement \[ \begin{align} \color{purple}{P(\theta_1,\theta_2~|~D)} =& \color{purple}{beta(\theta_1~|~z_1+a_1,N_1-z_1+b_1)} \\ & \color{purple}{beta(\theta_2~|~z_2+a_2,N_2-z_2+b_2)} \end{align} \]

On voit que lorsque le prior est le produit de deux distributions beta indépendantes, la distribution postérieure est le produit de deux distributions beta indépendantes

Représentation graphique

MCMC : modèle à deux variables - algorithme Metropolis

Calcul de la distribution postérieure par l'algorithme Metropolis

  • La proposition de saut repose sur une une loi normale 2D (généralisation ++)

\( ~\longrightarrow~ \) La position suivante peut-être n'importe où sur la grille

  • La règle d'acceptation est inchangée

MCMC : algorithme Gibbs - Description

Description de l'algorithme Gibbs

On voit que l'efficacité de l'algorithme Metropolis repose sur la proposition de déplacement.
Cette fonction peut-être difficile à calibrer:
\( \longrightarrow~ \) on voudrait que sd soit petit dans les zones à forte densité
\( \longrightarrow~ \) au contraire, on voudrait que sd soit grand dans les zones à faible densité

L'algorithme Gibbs résout ce problème

L'algorithme de Gibbs parcourt l'espace des paramètres comme l'algorithme Metropolis
Mais plutôt que de modifier tous les paramètres en même temps, il ne modifie qu'un paramètre après l'autre
\( \longrightarrow~ \) la proposition de déplacement est basé sur la fonction de densité conditionnelle
\( \longrightarrow~ \) De ce fait, le déplacement est toujours accepté

MCMC : algorithme Gibbs - Description

Algorithme Gibbs avec 2 paramètres


\( \longrightarrow~ \) Sélectionner un point de départ \( (\theta_{1init},\theta_{2init}) \)
On peut sélectionner n'importe quelle valeur

\( \longrightarrow~ \) Sélectionner le premier paramètre par exemple \( \theta_1 \)

\( \longrightarrow~ \) Générer une valeur aléatoire de \( \theta_1 \)
Faire un tirage directement suivant la probabilité \( P(\theta_1~|~\theta2,D) \)

\( \longrightarrow~ \) Générer une valeur aléatoire de \( \theta_2 \)
Faire un tirage suivant la probabilité \( P(\theta_2~|~\theta1,D) \)
Avec pour valeur de \( \theta_1 \) la valeur calculée précédemment

\( \longrightarrow~ \) On boucle sur les deux dernières étapes pour parcourir

illustration avec 2 paramètres

MCMC : algorithme Gibbs - Description

Calcul de la prochaine valeur du paramètre

Pour utiliser l'algorithme, il faut pouvoir calculer toutes les fonctions de densité conditionnelle

Dans notre exemple :

\[ \begin{align} P(\theta_1~|~\theta_2, D) &= \frac{P(\theta_1,\theta_2~|~D)}{\int P(\theta_1,\theta_2~|~D)~\mathrm d\theta_1} \\ &= \frac{beta(\theta_1~|~z_1+a_1,N_1-z_1+b_1)~beta(\theta_2~|~z_2+a_2,N_2-z_2+b_2)}{\int beta(\theta_1~|~z_1+a_1,N_1-z_1+b_1)~beta(\theta_2~|~z_2+a_2,N_2-z_2+b_2)~\mathrm d\theta_1} \\ &= \frac{beta(\theta_1~|~z_1+a_1,N_1-z_1+b_1)~beta(\theta_2~|~z_2+a_2,N_2-z_2+b_2)}{beta(\theta_2~|~z_2+a_2,N_2-z_2+b_2)~\color{orangered}{\int beta(\theta_1~|~z_1+a_1,N_1-z_1+b_1)~\mathrm d\theta_1}} \\ &= \frac{beta(\theta_2~|~z_2+a_2,N_2-z_2+b_2)}{beta(\theta_2~|~z_2+a_2,N_2-z_2+b_2)}~beta(\theta_1~|~z_1+a_1,N_1-z_1+b_1) \\ &= beta(\theta_1~|~z_1+a_1,N_1-z_1+b_1) \end{align} \]

En conclusion,
\[ P(\theta_1~|~\theta_2, D) = P(\theta_1~|~D) \]

MCMC : algorithme Gibbs - résultats


MCMC : algorithme Gibbs - Conclusion

  1. L'algorithme Gibbs est particulièrement utile lorsque la distribution jointe \( P(\{\theta_i\}~|~D) \) n'est pas calculable ou facile à échantillonner, mais que les densité conditionnelle \( P(\theta_i~|~\{\theta_{j \ne i}\}, D) \) le sont.

  2. L'algorithme Gibbs peut-être vu comme un cas spécial de l'algorithme Metropolis pour lequel la proposition de saut dépend de la position dans l'espace des paramètres et du paramètre sélectionné.

  3. En général, L'algorithme Gibbs améliore les performances de l'algorithme Metropolis

Planning

Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton

Plan

  • Rappels et objectifs
  • Introduction
  • Evaluation de la distribution posterieure par échantillonnage
    • L'algorithme Metropolis
    • L'algorithme Gibbs
    • L'algorithme Hamilton
  • Vérification du résultat

MCMC : algorithme HMC

Cet algorithme a été présenté pour la première fois en 1987
les auteurs :
Simon Duane, A.D. Kennedy, Brian Pendleton and Duncan Roweth





https://www.researchgate.net/profile/Simon_Duane

Méthodes Hamiltonian Monte Carlo

On a vu que l'algoritme Gibbs améliore les performances de convergence de l'algorithme Metropolis.
Mais il nécessite de pouvoir calculer les densités conditionnelles \( P(\theta_i~|~\theta_{j \ne i}, D) \)

L'algoritme est performant lorsqu'il y a indépendance des paramètres \( P(\theta_i~|~\theta_{j \ne i}, D) = P(\theta_i~|~D) \) Ce n'est pas le cas en général!!!

Pour le HMC, abandon de cette stratégie

Méthodes Hamiltonian Monte Carlo

L'algorithme Hamiltonian Monte Carlo est une amélioration de Metroplolis.

  • Une version particulière de l'algorithme Metropolis :
    \( \longrightarrow~ \) la proposition de déplacement est toujours une gaussienne multivariée centrée sur la position courante
    \( \longrightarrow~ \) la proposition de déplacement est une gaussienne qui garde toujours la même forme (variance constante)

  • L'adaptation de la proposition de déplacement en fonction de la position courante du paramètre.

Méthodes Hamiltonian Monte Carlo

Proposition de déplacement

\( \longrightarrow~ \) Sélectionner un point de départ
\( ~\bullet~ \) Il faut choisir \( ~\theta \in \left[-3,3\right] \) dans cet exemple


\( \longrightarrow~ \) Il faut définir une fonction de potentiel
On prend l'opposé du logarithme de la fonction de densité


\( \longrightarrow~ \) On peut générer un déplacement
On présente ici un ensemble de point pour voir la tendance


\( \longrightarrow~ \) L'histograme présente la tendance
On voit que le mode de la distribution de la future position n'est pas aligné avec la position du point de départ

Méthodes Hamiltonian Monte Carlo

L'algorithme HMC

\( \longrightarrow~ \) Sélectionner un point de départ
On peut sélectionner n'importe quelle valeur de \( \theta \) dans l'espace des paramètres

\( \longrightarrow~ \) Générer une nouvelle position à partir de la fonction de potentiel
\( ~\bullet~ \) L'algorithme génère une proposition par analogie au lancé d'une bille sur la fonction de potentiel
\( ~\bullet~ \) Le pas de déplacement, ainsi que le nombre de déplacement est fixé à l'avance
\( ~\bullet~ \) La force avec laquelle on lance la bille est déterminée par une loi normale centrée

\( \longrightarrow~ \) Accepter ou rejeter la proposition de déplacement
Faire un tirage suivant la probabilité suivante en considérant \( \phi \) comme le moment associé à la bille

\[ P_{déplacement}=min \left(\frac{P(\theta_{proposé}~|~D)~P(\phi_{proposé})}{P(\theta_{courant}~|~D)~P(\phi_{courant})},1 \right) \]

\( \longrightarrow~ \) on enregistre la nouvelle position et on recommence…

Méthodes Hamiltonian Monte Carlo


Faiblesse de l'algorithme


La faiblesse est liée à la génération de la proposition de déplacement.


quand le pas de déplacement est trop petit
ou
quand le nombre de pas est trop grand


Méthodes Hamiltonian Monte Carlo


Faiblesse de l'algorithme


La faiblesse est liée à la génération de la proposition de déplacement.


L'influence du premier potentiel choisi


Méthodes Hamiltonian Monte Carlo


  1. L'algorithme HMC est particulièrement utile lorsque le modèle contient un grand nombre de paramètres et quand les paramètres sont corrélés entre eux.

  2. Il nécessite un paramétrage assez fin. Des logiciels proposent des méthodes clés en main très efficaces.
    Exemple : STAN

  3. Le package brms (déjà vu aux chapitres précédents) utilise cet algorithme … \( \longrightarrow~~~ \) fonction : brm()

Planning


Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton


Plan

  • Rappels et objectifs
  • Introduction
  • Évaluation de la distribution postérieure par échantillonnage
  • Vérification du résultat
    • Représentativité
    • Précision
    • Implémentation

Vérifications du résultat des methodes MCMC

Pourquoi ces méthodes ne sont pas sûres à 100%

A. Le temps de calcul limité !

B. ces méthodes nécessitent un paramétrage fin
\( \longrightarrow~ \) Algorithme Gibbs
\( ~\bullet~ \) la variance de la distribution normale de la proposition

\( \longrightarrow~ \) Algoorithme HMC
C'est le paramétrage du calcul de la position suivante qui est délicat :
\( ~\bullet~ \) le pas
\( ~\bullet~ \) la longueur du pas
\( ~\bullet~ \) la force appliquée à la bille

Vérifications du résultat des methodes MCMC


Comment vérifier le résultat

Ces méthodes produisent des chaine de valeurs de paramètres (échantillons).
Chaque chaine produite doit respecter trois critères pour être validée

  1. Elle doit-être représentative de la distribution postérieure
    \( ~\bullet~ \) elle ne doivent pas dépendre du point de départ
    \( ~\bullet~ \) elle ne doit pas être cantonnée à une région de l'espace des paramètres

  2. Elle doit-être suffisamment longue pour assurer la précision et la stabilité du résultat
    \( ~\bullet~ \) La tendance centrale et le HDPI calculés à partir de la chaine ne doivent pas changer si on relance la procédure

  3. Elle doit-être générée de manière efficace

Vérifications du résultat des methodes MCMC


La représentativité

\( \longrightarrow~ \) Vérification visuelle des trajectoires
\( ~\bullet~ \) Les chaines doivent occuper le même espace
\( ~\bullet~ \) la convergence ne dépend pas du point de départ
\( ~\bullet~ \) Aucune chaine ne doit avoir de trajectoire particulière (cyclique)

\( \longrightarrow~ \) Vérification visuelle des densités
\( ~\bullet~ \) Les densités doivent se superposer
\( ~\bullet~ \) Les densités doivent-être régulières


Vérifications du résultat des methodes MCMC


La représentativité

Cette affichage ne montre que les 500 premières itérations

Les trajectoires ne se superposent pas au début (zone orange)

La densité est également affectée

C'est la Burn-in period

En pratique on supprime ces premières itérations


Vérifications du résultat des methodes MCMC


La représentativité

\( \longrightarrow~ \) Vérification numérique des chaines

le Shrink factor (en bas à gauche)
Correspond au “Rhat” de brms

C'est le rapport entre la variance inter chaines et intra chaines.
\( \bullet~ \) intuitivement la valeur doit tendre vers 1.0
\( \bullet~ \) en pratique elle est acceptable jusqu'à 1.1


Vérifications du résultat des methodes MCMC


Stabilité et précicion

Plus la chaine est longue et plus le résultat sera précis et stable
Si la chaine s'attarde sur chaque position, et que le nombre d'itrations reste le même, alors on perd en précision
Il lui faudra plus d'itérations pour arriver au même niveau de précision

\( \longrightarrow~ \) Mesure l'autocorrélation
L'autocorrelation est la corrélation de la chaine avec elle-même mais décalé de \( k \) itérations (lag)

\( \longrightarrow~ \) L'autocorrelation doit tendre vers 0 lorsque le lag augmente


Vérifications du résultat des methodes MCMC


Stabilité et précicion

La fonction d'autocorrélation est représentée pour chaque chaine (en haut à droite)

Un autre résultat rend compte de la précision de l'échantillon : l'effective sample size ESS

\[ ESS~=~\frac{N}{1+2\sum_k ACF(k)} \]

ESS représente la taille d'un échantillon non-autocorrélé extrait de la somme de toutes les chaines

Pour une précision raisonnable du HDI, il est recommandé d'avoir un ESS autour de 10000


Vérifications du résultat des methodes MCMC


Stabilité et précicion

L'erreur standard d'un échantillon est donné par le calul : \( SE~=~SD/ \sqrt{N} \)
plus N augmente plus l'erreur standard diminue

On peut généraliser cette idée aux chaines de markov

\[ MCSE~=~SD/ \sqrt{ESS} \]

Pour une précision raisonnable de la tendance centrale, il faut que cette valeur soit petite

Dans l'exemple, l'ESS est petite mais malgré ça le MCSE est petite et le mode est assez stable




Planning


Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton


Plan

  • Rappels et objectifs
  • Introduction
  • Evaluation de la distribution posterieure par échantillonnage
  • Vérification du résultat
  • Application

Application

L'implémentation d'un modèle douteux

set.seed(8)
y <- rnorm(100, mean = 0, sd = 1)
b5.1 <-
  brm(data = list(y           = y,
                  intercept_1 = 1,
                  intercept_2 = 1), 
      family = gaussian,
      y ~ 0 + intercept_1 + intercept_2,
      prior = c(prior(uniform(-1e10, 1e10), class = b),
                prior(cauchy(0, 1), class = sigma)),
      inits = list(list(intercept_1 = 0, intercept_2 = 0, sigma = 1),
                   list(intercept_1 = 0, intercept_2 = 0, sigma = 1)),
      iter = 4000, warmup = 1000, chains = 2,
      seed = 8)

Application

L'implémentation d'un modèle valide

set.seed(8)
y <- rnorm(100, mean = 0, sd = 1)
b5.2 <-
  brm(data = list(y           = y,
                  intercept_1 = 1,
                  intercept_2 = 1), 
      family = gaussian,
      y ~ 0 + intercept_1 + intercept_2,
      prior = c(prior(normal(0, 10), class = b),
                prior(cauchy(0, 1), class = sigma)),
      inits = list(list(intercept_1 = 0, intercept_2 = 0, sigma = 1),
                   list(intercept_1 = 0, intercept_2 = 0, sigma = 1)),
      iter = 4000, warmup = 1000, chains = 2,
      seed = 8)

Remplacement du prior uniform par un prior normal

Application

Affichage des chaines